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Nomenclature 


2t 


- wall  skin  friction  coefficient 


* * 

turbulent  eddy  viscosity  ■ 


* 

P 


- static  pressure 


* 

P 


<{> 


- body  surface  angle,  see  Figure  1 


- Stokes  stream  function  - 


P L U„ 


distance  normal  from  centerline  ■ — 


body  radius  ■ — j 
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Re 


Reynolds  number 


* * * 

P U L 
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velocity  component  In  x-dlrectlon 
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u 
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inviscid  velocity  at  body  surface 
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Nomenclature  (Continued) 


Reynolds  stress  in  boundary  layer 


* 

X 

X - arc  length  distance  along  body  ■ — j- 

L 

* 

y - distance  normal  to  body  surface  - ^ 

L 


Other  quantities  are  defined  in  the  text. 
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- reference  length 


- viscosity 


- kinematic  viscosity 


- density 


- wall  shear  stress  <■  y* 


3u 


Uy 


friction  velocity  ■ 


* 

IP  J 


U - free  stream  velocity 
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INTRODUCTION 

A second-order  accurate  Keller  box  solution  of  the  Incoapresslble 
axlsymmetrlc  boundary  layer  equations  has  been  enployed  by  Cebecl  [1] 
using  the  Mangler-Levy-Lees  transformation  to  avoid  a singularity  at 
the  stagnation  point.  Unfortunately  this  transformation  causes  a 
singularity  at  the  tall  of  the  body  making  extension  of  the  calculations 
Into  the  wake  region  Impossible. 

A second-order  accurate  code  has  been  written  by  Hoffman  [6]  which 
uses  non-transformed  variables  to  solve  the  equations.  This  code  Is 
Incapable  of  starting  at  the  stagnation  point;  therefore  a complete 
boundary  layer  profile  must  be  Input  from  the  Cebecl  program  mentioned 
above.  With  the  untransformed  variables  the  singularity  at  the  tall 
point  is  eliminated  and  better  accuracy  In  the  solution  Is  obtained 
near  the  tall. 

The  fourth-order  Keller  box  method  Is  a natural  extension  of  the 
second-order  sch'-se  presently  employed  by  Hoffman  (6) . The  term  fourth- 
order  as  used  here  refers  to  the  truncation  error  In  the  coordinate  normal 
to  the  axlsymmetrlc  body  surface.  The  truncation  error  In  the  tangential 
coordinate  is  still  second-order.  The  fourth-order  scheme  also  uses  the 
non-transformed  variables  mentioned  above  and  is  set  up  to  accept  a 
boundary  layer  profile  from  the  Cebecl  program  at  a station  downstream 
of  the  stagnation  point  In  the  same  way  as  the  method  of  Hoffman.  A 
fourth-order  solution  provides  higher  accuracy  everywhere  for  the  same 
mesh  spacing  or  provides  the  same  accuracy  for  a larger  mesh  spacing. 


t 

i 


I 

I 


I 
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when  compared  to  second-order  solutions.  Larger  mesh  spacing  requires  less 
computer  time  to  obtain  a solution;  therefore  the  fourth-order  solution  can 
provide  the  same  accuracy  with  less  computer  time  In  spite  of  the  Increased 
complexity  over  the  second-utder  method. 

INCOMPRESSIBLE  BOUNDARY  LAYER  EQUATIONS 

The  boundary  layer  equations  on  a body  of  revolution  for  steady  mean 
flow.  Including  transverse  curvature  effects,  are  [1]: 


where  dimensionless  variables  as  defined  In  the  nomenclature  have  been 
used  and  the  coordinate  system  Is  shown  In  Figure  1.  The  boundary 
conditions  are: 


u(x,  0)  - 0 , 

(3) 

V(x,  0)  - 0 , 

(4) 

11m  u(x,  y)  - u^(x)  , 

(5) 

y 00 

where  Ug(x)  Is  the  Invlscld  surface  velocity. 

It  Is  assumed  that  the  Reynolds  stress  Is  related 

to  the  velocity 

field  by  the  relation 

— i — r e 3u 

-u'v'  - — — 

Re  3y 

(6) 

C- 


1 


p 
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where  e is  the  eddy  viscosity.  An  algebraic  eddy  viscosity  model  is  used 
and  described  in  detail  in  Appendix  A.  Using  Equation  (6)  the  momentum 
equation  becomes 


— + ^ . 1 9 

9x  dx  r Re  9y 


r(l  + e) 

L J 


The  geometric  relation  between  r and  y from  Figure  1 is: 


(7) 


r(x,  y)  = TqCx)  + y cos  ip  . (8) 

The  term  y cos  p in  Equation  (8)  is  the  transverse  curvature  term  which 
is  important  at  the  aft  end  of  the  body  where  the  boundary  layer  thickness 
is  of  the  same  order  of  magnitude  as  the  body  radius. 

The  Reynolds  number  is  eliminated  from  Equation  (7),  in  the  case  of  a 
"thin"  boundary  layer,  by  introducing  the  scaled  variables: 


I 

f 


i 


I 

1 


Re 


Re 


Equations  (1),  (7),  and  (8)  become  respectively. 


(ru)  + (rv)  - 0 , 

9y 


9u  — 9u  dp.  19  F/il  \9u~I 

+ V—  = ^ — r(l  + e)  — 

3y  ' 3y  L 


r(x,  y)  - r-(x)  + — ^ cos  i}) 

/Te 


(9) 


(10) 


(11) 


(12) 


(13) 
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In  order  to  use  the  box  method,  these  equations  must  be  rewritten  In 
first-order  form.  The  Stokes  stream  function  Is  Introduced,  defined  by 


ru  ■ — ^ 

3y 


- -Dll) 

■V  • ■* 


rv 


Dx 


(14) 


The  continuity  equation.  Equation  (11),  la  Identically  satisfied  by  the 
stream  function.  The  Euler  eqautlon  at  the  body  surface  which  relates 
the  pressure  gradient  to  the  Invlscld  surface  velocity  Is: 


, du 

dp  e 

- -u  -3 — 
dx  e dx 

If  use  Is  made  of  the  laminar  shear  stress,  defined  by 


3u 

3y 


(15) 


(16) 


then  the  momentum  equation.  Equation  (12),  can  be  written  In  the  following 
first-order  form: 


r 3 . 2 2,  ^ ^ 3 f..  _ - 

3y 

where 

b - 1 + e . 

Summarizing,  the  set  of  first-order  equations  to  be  solved  Is: 


r3.2  2.  ^ ^ 

7=  ” 2 ^ - “e  ^ ^ ' 

3y 


W 

3y 


^(U) 


ru 


(17) 


(18) 


(19) 


(20) 


(21) 


3y 


1 

1 
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The  new  boundary  conditions  become  [6], 

u ■ “ 0 at 

y - 0 , 

(22) 

u u as 

e 

y -►  00  . 

(23) 

DEVELOPMENT  OF  FOURTH-ORDER  EQUATIONS 

The  finite-difference  procedure  used  in  solving  the  parabolic  boundary 
layer  equations  is  the  box  method  of  H.  B.  Keller  [2].  The  first-order 
form  of  the  governing  differential  equations,  (19)-(21),  when  put  in 
finite  difference  form  will  result  in  algebraic  equations  which  relate 
unknowns  along  line  (i+1) , as  shown  in  Figure  2,  to  line  i . The  solution 
along  line  1 is  assumed  known  and  since  the  equations  are  parabolic  the 
solution  at  line  (i+1)  can  be  obtained.  The  result  is  a marching  procedure 
beginning  at  the  initial  station  and  proceeding  downstream. 

In  applying  Keller's  box  method  one  approximates  the  equations  by 
difference  quotients  about  point  P , the  geometric  center  of  the  box 
element,  as  shown  in  Figure  2.  The  derivatives  in  the  axial  direction 
are  second-order  accurate  and  are  given  by  simple  centered  differences. 

The  fourth-order  accurate  expression  used  for  derivatives  in  the  normal 
coordinate  is  given  by  Wornom  [3]  as 


2 

*j+l  ~ ~ ~2  [®j+l  (®j+l"  " ®j"]  ^ ^ ° ’ 


where  primes  denote  differentiation  in  the  normal  direction  (w.r.t.  y) 
and  hj  is  the  step-size  in  the  normal  direction  given  by 
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g represents  a matrix  of  the  differentiable  quantities  from  Equations 
(19)-(21)  given  by 

(26) 


r>  . 


rbx 


u 


Ftom  these  equations  the  first  derivative  of  g is  obviously 


^ . r 3 .2 

3^  ^ 2 ^ - 


“e'> 


ru 

T 


(27) 


The  form  of  g"  , however,  is  not  so  obvious,  and  becomes  quite  Involved 
algebraically.  In  differentiating  Equation  (27)  one  obtains  terms  which 
contain  the  derivative  t'  . Here  the  turbulence  model  used  to  determine 
the  eddy  viscosity  e becomes  of  great  significance.  In  the  inner  region 
e is  given  by 

e - 1t1  , (28) 

where 

C - Re^^^  — , (29) 

and  where  £ is  the  mixing  length,  also  a function  of  y . In  the  outer 
region,  however,  e is  constant  and  independent  of  boty  y and  T . See 
Appendix  A for  a more  detailed  description  of  the  turbulence  model. 

The  derivative  T*  may  be  determined  from  Equation  (19)  by  using  the 
product  rule  of  differentiation  and  Isolating  T*  . In  the  outer  region 
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the  value  of  b from  Equation  (18)  Is  constant  and  is  left  as  such.  However, 
in  the  inner  region  the  value  of  b is  a function  of  y , and  must  also  be 
differentiated.  Thus  the  expression  for  t'  in  the  inner  region  becomes  much 
more  complicated  than  that  in  the  outer  region.  Appendix  B shows  the  deriva- 
tion of  g"  for  both  of  these  cases,  the  results  of  which  appear  below: 


- 5 ^ ^ ^ . a 3 . 2 2-  1 1 

_x  _ (ru)  + (ru)  — + (u  - u^  ) - 

n.  3 , 2 2.  T ^ .2  a , a ^fl 

hr7^(u  -u  ) ■::^-s't T ST 

[2  3x  e r3x  r 


g inner 


{ au  + rt  } 


f 1 1 n.  3 . 2 2,  T 3^  .2  a ^ a ^2] 


( 

f '\ 

^ ^ \ ^ / v3T.a3,2  2.  ^l.r 

-T  ^ (ru)  + (ru)  — + - — (u  - u^  ) - j 


g outer  = 


9/2  2, 

(u  - u ) - abx  - T 
3x  e 


{ au  + rt  } 

1 [r  3 / 2 2.  . 3iin 

7^  [I  37  - “e  > - ^ • 


Equations  (26),  (27),  and  either  (30)  or  (31)  are  then  substituted  into  the 

finite-difference  expression.  Equation  (24)  to  obtain  three  nonlinear  finite- 

0 / \ 

difference  equations.  The  — r — terms  in  these  equations  are  approximated  by 
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a central  difference  quotient,  as  already  mentioned,  making  the  equations 
second-order  accurate  In  the  x coordinate  and  fourth-order  accurate  In 
the  y coordinate.  The  non-linear  teims  are  then  linearized  by  Newton's 
method  to  obtain  three  linear  finite-difference  equations, developed  In 
Appendix  C,  which  are  shown  below: 


AAA 

^^i+i,j^’*'i+i,j  ®^i+i,j^“i+i,j  ^^1+1, 1+1, j 


A 


A 


A 


°^i+l,j'^'^i+l,j+l  ^^i+l,j'^'^i+l,j+i*‘  ^^i+l,j'^^i+l,j+l“  ^^1+1, j ’ 


^^i+l,j^’^i+l,j  ®^i+l,j^i+l,j  ^^1+1, 1+1, j °^i+l,j^'^l+l,j+l 


A 


A 


A 


®^^i+l,j*^“i+l,j+l'’'^^i+l,j'^’^i+l,j+l  ^^1+1, j ’ 


(33) 


A 

^^i+i,j‘^“i+i,j+i 


A A 

^^i+l,j^'^i+l,j+l  ~ ^^1+1, j 


(34) 


where  the  "6"  terms  in  the  above  equations  indicate  corrections  to  , 
u , and  T . The  appropriate  boundary  conditions  become 


^n+1,1 

"^“i+l,! 

*^“i+l,N+l 


0 

0 

0 


(35) 

(36) 

(37) 


1 ' ...  -i.- 
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MATRIX  FORMULATION  AND  SOLUTION 

In  order  to  put  the  equations  in  trldlagonal  matrix  form,  the  boundary 
Conditions  given  by  Equations  (35) -(37)  and  the  finite-difference  equations 
given  by  Equations  (32) -(34)  are  grouped  in  the  following  fashion  (i+1  sub- 
scripts are  understood) : 


= 0 


6u  = 0 

A A 


A 

A 

A 

A 

A 

1 

+ C3j^6t^  + D3j^6t2 

+ £3^^602 

+ F33^6t2 

• s^i 

j 

1 

A 

A 

A 

A 

A 

i 

1 

+ + D1j^«S4>2 

+ E1j^6u2 

+ F12^5t2 

1 

1 

A 

A 

A 

A 

A 

. 2*'*^  ! 

^ 1 

+ C2j^6Tj^  + D2j^6ijj2 

+ E23^6u2 

+ F23^6t2 

= S2^ 

A 

A 

A 

A 

A 

+ C326T2  + D326<I'3 

+ E326u^ 

^ ^^26X3 

-52^ 

i 

A A A A A A A 

A^j-l^^-l  ®"j-l«“j-l  = sij. 

A , A A A A A A 

+ C2._^6Tj_^+  D2j_^6.Pj  + E2j_^6u^  + F2._^6Tj  = S2., 

A ^ A A A A A A 

+ B3j6uj  + C3j6Tj  + 03^6^^^^  + E3j6Uj^^  + f3j6Tj^^  - S3. 


A 

A 

A 

A * 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

+ B2j^% 

^^n'^^N+1  ‘ 

= S2, 

(N+1) 


^“n+1  “ 0 
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(38) 

the  preceding  equations  can  be  written  as  the  following  system  of  matrix 
.nations: 

Vi  ^1^2  “ h • 

®N+l^N  \+l^N+l  “ ^+1  ’ 


By  defining  a 3-component  column  vector  Zj  according  to 
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A n 

I 

i 

1 t 

f 

^+1  “ 

A A 

A 

F2 

N 

\ 

, (46)  i 

0 1 

0 

'( 

t 

i 

! 

and  the  are  3-componenc  column  vectors  as  follows: 


0 

0 

S3, 


(47) 


A ~1 

^^J-1 


A 

S2 


j-1 


A 

S3, 


(48) 


and 


^1, 


A 

S2, 


(49) 


Then  Equations  (39)-(41)  can  be  written  In  the  compact  block  matrix  form: 

Aj.  Z * R , (50) 

where  is  an  (N+1)  x (N+1)  block  trldlagonal  matrix  given  by 


h Cl 


^2  ^2  ^2 


C 


'n.\\ 

®N  N 


"N+1 


\+l 


(51) 


! i 


I 
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and  Z and  R are  (N+1) -component  column  vectors  as  follows 


Because  the  coefficient  matrix  is  block  tridiagonal.  Equation  (50) 

can  be  solved  by  an  efficient  elimination  procedure  similar  to  that 
used  by  Keller  [10],  This  procedure  is  presented  in  detail  in  Appendix  D. 
The  result  is  a solution  for  each  of  the  correction  terms  6u,  and  6t 
in  Equations  (32) -(34). 

These  correction  terms  are  then  added  to  the  values  of  , 

1+1  >J 

u^^j^  j+1*  order  to  obtain  new  iterates  for  these  variables. 

The  iteration  procedure  is  continued  until  prescribed  convergence  toler- 
ances are  met. 


RESULTS 

In  this  paper  computer  solutions  using  the  present  fourth-order  box 
method  are  compared  to  the  second-order  box  solution  of  Hoffman  [6]  and 
with  the  experimental  results  of  Patel  and  Lee  [14]  for  a low  drag  body 
revolution.  The  body  used  was  the  F-57  body,  shown  in  Figure  3;  its 


I 
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coordinates  are  given  below  [14] ; 


For  0 < X-  < X (fore-body) 
— u — m 


m 


-1.1723  + 0.7088  5^  + 1.0993  + 0.3642 


(53) 


For  x^  ^ Xq  1 (pointed  aft-body)  , 


-0.11996  ^2  - 2.58278  + 3.52544  + 0.17730  [ . (54) 


where  ^ = x^/x  , = (1  - x„)/(l  - x ) , x„  is  the  axial  distance  measured 

1 U m z 0 m 0 

from  the  nose,  r is  the  local  radius,  x is  the  axial  location  of  the  maxi- 
o m 

mum  radius  r , and  all  length  variables  have  been  non-dlmensionalized  with 
m 

respect  to  L , the  body  length  (L  = 1.2l9m).  Thus,  x = x /L  * 0.4446, 

m m 

* / * 6 
r = r /L  = 0.1170.  All  solutions  were  at  a Reynolds  number  of  1.2  x 10 

” ” U*  L* 

00 

where  Re  = 


. The  experimental  measurements  [14]  were  taken  with 


* 1 * 

the  boundary  layer  tripped  at  Xq  = Xq  /L  = 0.475  by  a circular  trip  wire 
of  1.664  mm  diameter  wrapped  around  the  body.  Subsequent  analysis  of  data 
revealed  that  the  downstream  influence  of  such  a relatively  large  trip 
wire  may  have  persisted  up  to  Xq  - 0.6  [14] . The  turbulence  model  used 
in  the  computer  code  cannot  accurately  simulate  such  a sudden  jump  from 
laminar  to  turbulent  flow;  thus  a transition  Intermlttency  factor  has 
been  Included  as  described  in  Appendix  A.  The  axial  distance  where  the 
transition  region  should  begin  in  the  computer  solution  was  then  chosen 


r 
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by  equating  the  momentum  deficit  area  of  the  calculation  to  that  of  the 
experiment.  Following  Hoffman  [11]  an  axial  station  at  x^  ■ 0.3271  was 
used  In  the  calculations. 

The  pressure  distribution  Input  to  the  computer  programs  was  taken 
from  the  strong-interaction  calculations  of  Hoffman  [11].  The  Iterative 
procedure  used  to  generate  these  data  Is  based  on  a correction  of  the 
displacement  body  Idea  to  represent  streamline  curvature  effects  by 
means  of  a simple  pressure  mapping,  and  Is  described  In  Reference  11. 
Figure  A shows  this  distribution  as  a function  of  the  axial  body 
coordinate. 

As  mentioned  previously,  in  order  to  avoid  a singularity  at  the  nose 
point,  the  computer  codes  for  both  the  second-  and  fourth-order  box 
solutions,  which  use  non-transformed  variables,  are  set  up  to  accept  as 
input  a boundary  layer  profile  at  some  station  downstream  of  the  nose. 

The  station  was  chosen  in  the  laminar  region  at  Xq  = 0.09548  and  the 
profile  was  generated  with  Hoffman's  second -order  computer  solution 
which  uses  the  Mangier -Levy-Lees  transformation  [6] . 


Integral  Parameters 

Calculations  were  made  of  the  integral  parameters  C^,  0,  and  H at 
each  axial  coordinate  Xq  where  0 is  the  momentum  flux  deficit. 


00 


J 

0 


(55) 
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H is  the  boundary  layer  shape  factor. 


6 Is  the  mass  flux  deficit 


00 


O 


and  is  the  wall  friction  coefficient. 


(56) 


(57) 


\ . (58) 

P 

Results  are  presented  in  graphical  form  in  Figures  5-7  for  C^,  0,  and 
H, respectively,  as  functions  of  x^.  The  fourth-order  calculations, 
which  start  with  an  input  boundary  layer  profile  of  11  points  in  the 
normal  coordinate  at  an  axial  station  of  Xq  • 0.09548,  are  compared 
with  the  second-order  calculations  [6],  starting  with  an  input  boundary 
layer  profile  of  31  points  in  the  normal  coordinate  at  the  same  axial 
station.  The  experimental  results  [14]  are  also  shown  for  comparison. 
Figure  5 shows  that  the  calculated  agrees  well  with  experiment  for 

0.6  < Xq  < 0.8,  but  disagrees  sharply  both  fore  and  aft  of  this  interval. 
From  Figure  6 one  can  see  that  the  calculated  0 is  lower  than  experi- 
ment for  the  interval  0.6  < x^  < 0.9,  but  slightly  higher  aft  of  this 
interval.  The  calculated  shape  factor  H in  Figure  7 has  very  poor 
agreement  with  experiment.  The  reason  for  these  disagreements  is 


I 
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attributed  to  the  inability  of  the  algebraic  turbulence  model  to  accurately 
simulate  the  large  trip  wire  used  in  the  experiment,  and  to  the  general 
Inadequacy  of  the  turbulence  model  In  the  tall  region. 

All  these  figures  show  that  agreement  between  second-  and  fourth- 
order  calculations  is  quite  good;  the  differences  between  the  two  being 
Indistinguishable  on  most  of  the  body. 

Boundary  Layer  Profiles 

Figures  8-10  show  the  boundary  layer  profiles  at  axial  coordinates 
of  0.601,  0.880,  and  0.990, respectively.  The  calculated  profile  at  x^  ■ 0.601 
(see  Figure  8)  reflects  the  Influence  of  the  large  trip  wire  used  in  the 
experimental  claculatlons;  the  calculated  results  near  the  wall  are  less 
full  than  the  experimental  results.  At  the  two  downstream  stations  where 
the  effect  of  the  large  trip  wire  has  died  out,  the  opposite  Is  true;  the 
profiles  are  slightly  too  full  near  the  wall,  as  can  be  seen  In  Figures  9 
and  10,  for  x^  = 0.880  and  0.990  respectively.  These  discrepancies  are 
attributed  to  the  Inadequacy  of  the  algebraic  turbulence  model  near  the 
tail  of  the  body. 

In  Figures  8-10  the  second-  and  fourth-order  calculations  are  seen  to 
be  Indistinguishable.  Note  that  the  fourth-order  calculations  were  made 
with  a mesh-spaclng  three  times  as  large  as  that  of  the  second-order 
calculations . 

Computer  Run  Times 

As  mentioned  previously,  both  the  second-  and  fourth-order  calculations 
were  performed  on  the  Penn  State  IBM  370/3033  computer  system.  The  net 
C.P.U.  time  for  the  second-order  calculations  shown  In  the  figures  was 
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38  seconds,  whereas  the  net  C.P.U.  time  for  the  fourth-order  calculations 
with  one-third  of  the  number  of  points  in  the  input  profile  was  28  seconds. 
This  represents  a 26  percent  decrease  in  run  time  for  the  fourth-order 
method. 

CONCLUDING  REMARKS 

A fourth-order  box  solution  of  the  incompressible,  axisymmetric, 
laminar  or  turbulent  boundary  layer  equations  has  been  presented  and, 
as  an  example,  applied  to  a low-drag  body  of  revolution.  For  this  case, 
the  fourth-order  solution  had  an  input  mesh-spacing  three  tljoes  larger 
than  that  of  a second-order  solution,  but  provided  results  of  equivalent 
accuracy  and  with  26  percent  less  computer  run  time.  Fourth-order  box 
solutions  of  other  equations  have  also  been  found  to  provide  the  equiv- 
alent accuracy  of  second-order  solutions  with  fewer  mesh  points  and 
requiring  less  computer  time.  For  example,  Wornom  [3]  has  achieved  such 
a result  for  the  two-dimensional  stagnation  point  solution. 

One  disadvantage  of  the  fourth-order  box  solution  presented  here 
should  be  noted.  The  solution  sometimes  has  greater  convergence 
difficulties  at  the  outer  edge  of  the  boundary  layer  than  does  the 
second-order  solution  with  the  same  input  conditions.  Convergence 
problems  have  been  noted  for  solutions  with  very  small  mesh-spacing 
or  with  smaller  convergence  tolerances.  These  problems  were  especially 
prevalent  when  the  one  piece  turbulence  model  of  Clowacki  [12]  was  tried. 

No  convergence  problems  were  encountered  with  the  second-order  method 
using  this  model,  but  the  fourth-order  method  would  not  converge  to  a 
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solution  aft  of  axial  coordinate  Xq  ••  0.880.  The  reason  for  this  behavior 
Is  not  knovm;  one  possibility  Is  that  because  of  the  Increased  complexity 
of  the  fourth-order  coefficients,  numerical  noise  Is  amplified  as  these 


coefficients  get  very  small  at  the  outer  edge  of  the  boundary  layer. 

This  problem  requires  further  Investigation  . 

Other  fourth-order  accurate  methods  are  available  which  have  not 
been  considered  In  this  paper.  See,  for  example,  a comparison  of  higher- 
order  solutions  by  Rubin  and  Khosla  [13].  On  the  other  hand,  spline 
methods  may  provide  the  same  advantages  as  the  fourth-order  Keller  box 
method,  but  with  much  less  algebraic  manipulation  (a  quantity  of  great 
abundance  In  the  method  presented  In  this  paper) . 
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Appendix  A 

Algebraic  Eddy  Viscosity  Model 
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This  appendix  describes  the  details  of  the  algebraic  eddy  viscosity 

model  and  treatment  of  the  transition  region.  In  this  model  the  boundary 

layer  is  divided  into  an  inner  (law  of  the  wall)  and  an  outer  (law  of  the 

wake)  region;  each  region  has  its  own  equation  with  the  junction  between 

the  two  occurring  when  e.  = e . 

1 o 

Inner  Region 

The  inner  eddy  viscosity  used  is  that  of  Crawford  and  Kays  (5]  as 
modified  by  Hoffman  [6]  for  a thick  axisymraetric  boundary  layer.  The 


final  form  is 


* 'o  97 

where  £.  is  the  mixing  length  given  by 


(A.l) 


i = 0.41  Y 


[l  - exp  (-  I)]  , 


(A. 2) 


and  Y is  the  modification  for  thick  axisymmetric  boundary  layers  mentioned 
above  and  is  defined  as 


Y ^ In  i- 

• cos*  r 

o 

A is  the  sublayer  thickness  parameter;  in  universal  (law  of  the  wall) 


(A.3) 


variables, 


A " A Re  u 


(A.4) 


where  u^  is  the  dimensionless  friction  velocity  given  by 


“t  “ jCf 


(A.  5) 


and  is  the  skin  friction  coefficient  based  on  U^.  The  empirical 
00 

formula  of  Crawford  and  Kays  [5]  is  used  for  the  sublayer  thickness 


parameter,  which  for  an  Impervious  wall  is 
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^ 1.0+  aP 

where  A is  the  flat  plate  value  of  = 26.0, 

[30.175  if  < 0 


® ] 20.590  if  P"^  > 0 , 


u du 
e e 


(A.  6) 


(A.7) 


(A. 8) 


4. 

The  above  expression  for  A was  derived  from  data  for  equilibrium 
two-dimensional  boundary  layers.  In  a rapidly  changing  pressure  gradient 
the  boundary  layer  does  not  Instantaneously  adjust  to  a new  equilibrium 
state.  Instead,  it  tends  to  lag  in  its  adjustment;  to  allow  for  this  lag 
Crawford  and  Kays  [5]  have  Introduced  a lag  equation  given  by 


p-^-p-^ 


+ 4* 

where  T is  a time  constant,  P is  the  effective  value  of  P to  be 

ett 

used  in  the  relation  for  A (Equation  A. 6),  and  P is  the  local  value 
of  P^  as  given  by  Equation  (A. 8).  The  value  of  T used  is  Hoffman's 
modification  [6]  of  Crawford  and  Kay's  suggestion  [5]  which  is 


(A.9) 


(A. 10) 


Equation  (22)  is  solved  by  using  centered  differences  and  defining 

T 

to  get  the  finite-difference  form 


1 ff  + “+  ‘ 

1 [1^  ' 2tJ  ^ eff^  ^i+1/2 


(A. 11) 


(A.12) 


Equation  (A.12)  is  an  iterative  equation  because  of  the  term  so  that  each 

time  the  vector  of  unknowns  at  (i+1)  is  updated  in  solving  the  non-linear 
finite-difference  equations  of  motion.  Equation  (A.12)  must  also  be  updated. 
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Outer  Region 


In  the  outer  region  of  the  boundary  layer  the  eddy  viscosity  is 
taken  to  be  constant  using  the  form  suggested  by  Clauser  which  when 
non-dimensionalized  becomes  [7] 


e = Re  a u 6, 
o e k 


(A.13) 


where  is  the  kinematic  displacement  thickness  given  by 


(A. 14) 


and  o is  a dimensionless  quantity  which  accounts  for  low  Reynolds  number 


effects.  Cebicl  has  obtained  the  following  expression  for  a [ 8] : 

1 + IT 

o 

o » a -r  ■ — 

O 1 + TT 


where 


a - 0.0168 
o 


TT  - 0.55 
o 


TT  * TT^  [1  - exp  (-0.243  /~Z“  -0.298Z^)]  , 

Re„ 


(A.15) 


(A.16) 


(A.17) 


(A. 18) 


(A. 19) 


and  Re^  is  the  kinematic  momentum  thickness  Reynolds  number  defined  by 


Reg  - Re 
k 


(A. 20) 


where 


0 - — f — fl  - — 1 dy 

^ I,  “e  ^ 


(A.21) 


Note  that  when  Re.  < 425,  the  argument  of  x in  Equation  (A. 18)  becomes 
®k 

imaginary.  For  this  case  one  sets  ir  to  zero. 


) 

f 
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Transition  Region 

Since  an  algebraic  eddy  viscosity  is  used  to  model  turbulence  effects 
in  the  boundary  layer,  there  is  no  built-in  transition  region.  From  an 
empirical  curve  fit  discussed  in  Reference  6 a point  is  picked  along  the 
body  where  turbulence  is  "switched  on,"  or  tripped.  It  is  not  realistic 
to  model  the  transition  region  by  simply  preceding  from  laminar  to 
turbulent  flow  in  one  streamwise  step.  Therefore  to  account  for  the  finite 
distance  required  in  the  transition  region  for  the  boundary  layer  to  become 
fully  turbulent,  an  intermlttancy  factor  is  used  to  reduce  the  intensity 
of  the  turbulence  parameters.  This  factor,  from  Chen  and  Thyson  [9],  is 
given  by 

Ytr  ■ 1 - exp(-g  u^^  I2)  , (A.22) 


where 


Re 

1200  X. 


-1.34 


tr 


X 


(A. 23) 

(A. 24) 


(A. 25) 


Re  = Re  u x^  . (A. 26) 

Xtr  e tr 

The  trapezoidal  rule  is  used  in  the  computer  code  to  evaluate  the  running 
Integrals  and  1^.  The  eddy  viscosities  in  the  transition  region  are  then 
miltiplied  by  y^^  to  give 


Appendix  B 

Differentiation  of  Vector  g' 
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This  appendix  describes  in  detail  the  differentiation  of  the  vector  g'  with 
respect  to  y to  find  g"  in  both  the  inner  and  outer  regions  of  the  boundary  layer. 
Inner  Region 


Equation  (27)  can  be  written  in  component  form  as 

I j.  r 3 / 2 2. 


(B.l) 


g2  “ ru 


Differentiating  Equation  (B.l)  yields 


(B.2) 

(B.3) 


e"  = ' T - i r'  ^ . r [3  . 2 2. 

h laxj  ^ 3x  ^ ^ 2 3x  ■ “e  ^ ^ 2 "'^e  ^ J . (BJ 

Using  Equations  (20)-(21)  simplifies  the  above  equation  to 

*l"  + . ■ (B.5) 


where  from  Equation  (13) , 


a » r’  - ^ 


The  second  component  of  g"  becomes 
82"  “ au  + rx  , 
and  the  third  coinponent  becomes 


g,"  = t' 


To  evaluate  x'  the  product  rule  is  applied  to  Equation  (19)  to  get 

rbx'  + rxb’  + obx  = ^ (u^  - u ^)  - t 

2 3x  ' e 3x 

Differentiating  Equation  (18)  and  substituting  Equation  (28)  yields 


where 


b'  =■  s'x  + sx' 


3 “ t/|x| 


(B.6) 


(B.7) 


(B.8) 


(B.9) 


(B.IO) 


(B.ll) 


J 
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Thus,  by  substituting  Equation  (B.IO)  into  (B.9)  and  isolating  x', 
the  following  expression  is  obtained: 


e " = x'  = 

^3  (1+2 


1 1 & . 2 2.  X 2 a cx  2 

— ) ® T X sx  (B.12) 

sx)  23x  er3x  r r 


The  value  of  x obtained  above  must  be  substituted  into  Equation  (B.5)  to 


,1  3 f \ ■ / \ 3x  , a 3 / 2 2.  3il» 

_1  1 fl  3 , 2 2.  X 3i|)  ,2a  a 2 

l+2sxj  [2  3x  “ “e  ^ ■ r 3x  ■ ® " r “ r J * 

Combination  of  Equations  (B.13),  (B.7),  and  (B.12)  yields  the  final 

form  of  ss  given  by  Equation  (30)  . 

Outer  Region 

In  the  differentiation  of  g'  in  the  outer  region  it  becomes  apparent 
that  Equations  (B.5)  and  (B.7)  still  hold.  A different  expression  for 
t'  must  be  developed  however  because  the  eddy  viscosity  in  the  outer  region 
is  now  a constant.  From  Equation  (18)  it  is  srsen  that  b is  also  a constant; 
Equation  (19)  then  becomes 

rbx'  + obx  = 7 (u^  - u ^)  - x , 

2 3x  e 3x  ’ 


which  leads  to 


" , 1 fr  3 , 2 2.  3t{» 

®3  ’ ^ “ rt  [2  3^^  - “e  ^ 37  • 

Substitution  of  Equation  (B.14)  into  (B.5)  results  in  the  following 

expression: 

II  3 f N j.  / \ ^x  , a 3 . 2 2.  dilj  fll  fr  3 

h ■ Sr  <'“>  + 2 <"  - “e  > - sJ  [rtj  [2 

(u^  - u ^)  - abT  - T 

e 3x 


(B.14) 


(B.15) 


Combination  of  Equations  (B.15),  (B.7)  and  (B.14)  yields  the  final  form  of 


t 

1 I 


^ "outer  given  by  Equation  (31). 
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Appendix  C 

Derivation  of  Finite-Difference  Expressions 


This  appendix  describes  in  detail  the  derivation  of  the  three  linear 
finite-difference  expressions  given  by  Equations  (32)-(34). 

The  first  step  in  the  derivation  is  to  substitute  Equations  (26),  (27), 
and  either  (30)  or  (31)  into  Equation  (24)  to  obtain  a set  of  first-order, 
nonlinear  finite-difference  expressions  for  either  the  inner  or  outer 
regions  of  the  boundary  layer.  The  inner  layer  is  considered  first;  these 


expressions  become 

It  - “e^>  -'j  h (f)j  h '“j'  - %'>]  * if  h 
- “e')  - (7)^,2  ■ t?  - 

('  "Ij+J  " ^ ^ ‘'3’  - If  ^ 

h <»i>  [if  ("5'  - “ei-  (7)2  f - [7  - 

; ST^  \ - 0 , (C.l) 


"^j+l  9x  ^’*'3+1^  2 j+1 


j+1  3x 


2 2. 

^j  - “e  > + 


(C.l) 


’"j+l  ■ *^3  " 


jj  ^ 12  [' 


auj^l  + (rT)j^^ 


‘1  - 1 


(C.2) 


^ I 


I 
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and 


Consider  Equation  (C.2)  first.  Since  there  are  no  x-derivatives  it  can  be 
evaluated  at  (i+l,  j+l/2)  which  corresponds  to  point  Q in  Figure  2.  Equations 
(C.l)  and  (C.3),  however,  contain  x-derivatives  which  are  approximated  to 
second  order  by  using  central  differences,  and  the  expression' is  evaluated  at 
point  P in  Figure  2,  (i+1/2,  j+l/2),  the  center  of  the  box  element.  The 
resulting  expressions  are  quite  lengthy  and  will  not  be  shown  here.  However,  some 
examples  of  terms  which  appear  in  Equations  (C.l)  and  (C.3)  appear  below: 


3x  Ax^'  ^'•'l+l.j  ■ '''i,j^’ 

ot  2 1 f.  I 2.  / * 2. 

^ i+1/2, j+1  2 ^i+l,j+l+^®  ^i,j+l 


(C.4) 


(C.5) 


As  stated  previously,  the  boundary  layer  equations  are  parabolic, 
that  is,  there  is  no  upstream  influence  so  that  the  solution  is  obtained 
by  marching  downstream.  Thus,  only  the  values  along  the  station  (i+l) 
are  unknovm.  The  above  expressions  are  nonlinear  and  must  be  linearized 
in  order  to  obtain  a solution.  Since  Newton's  method  is  used  to  solve  the 
finite-difference  equations,  the  expressions  are  linearized  about  known  values 
at  station  i and  the  known  boundary  conditions  at  station  (i+l),  viz. 
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(1)  - ... 


•^i+l.N+l 


2 S j < N , 


(C.6) 


'i,N+l  ’ 


1+1,1  ^ 

(1) 

'l+l,j  " “i.j  ’ 

*1+1, N+1 


2 < j < N , 


(C.7) 


1+1,1 

'l+l,j 


T.  . , 2 < j < N 

l.J  i 


(C.8) 


T - T 

1+1, N+1  ’^1,N+1  ’ 


where  j = N+1  represents  the  grid  point  at  the  outer  edge  of  the  boundary 
layer  and  the  superscript  denotes  the  Iteration  count.  The  higher-order 
Iterates  are  then 


4+1,  j 

( 

4+1,  j 

I 

4+1,  j 


(n+1) 
' (n+1) 
' (n+1) 


4+1,  j 

I 

4+1,  j 

( 

4+1,  j 


4+1,  j 

I 

4+1,  j 

I 

4+1,  j 


, ) n ••  1,  2,  ...  (C.9) 


which  arc  Introduced  into  the  finite-difference  equations  developed  from 
Equations  (C.1)-(C.3),  with  only  linear  terms  in  the  corrections  being 
retained  In  the  result.  After  much  algebra,  which  will  not  be  shown 
here,  a coupled  linear  system  for  the  corrections 

obtained  for  1 i j i N+1,  which  is  represented  by  Equations 
(32)-(34).  The  superscript  n is  understood  in  these  expressions,  and  the 


coefficients  are  given  by 


I 
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Al. 

= -2 

i-H.j 

*■  > 

^•+1 
1+1,  j 

* ’^l+l.J  ■ ~6~  ^“i+1^  ’ 

A 

Cl  , 

±+1,5 

6 ’'1+1,  j 

A 

= 2 

±+1,5 

f 

^1+1.3 

“ "*'j  ’^i+l,j+l  “6“  ^“i+1^  * 

2 

A 

FI,  , 

SS  ' *■  1 ■ 

i+l,j 

6 ^i+l,j+l  ’ 

(C.]0) 

(C.ll) 

(C.12) 

(C.13) 

(C.14) 

(C.15) 


A 

SI 


1+1  ^“±+l.j+l  “l+l,j^  ^"^i+l.j+1  “ ^’'^^i+1 


.jj  ’ 


(C.16) 


A 

A2 


(t. 


i+l,j  Ax  '‘i+l,j  ‘i.j 


+ T.  .)  + 


h.‘ 


^ 1 ' 


Ax 


+ (Ha) 


i+1 
2 


..]■ 


(C.17) 


(C.18) 


i+l.j  “ 4-  [^““^i+l.j  [^]  - ^"^>i+l.j 

‘^^i+l.j  “ -2^+1^  ^®">l+l,j  ('^i+l.j  - -"-6“  (zr) 

1 X 

[-<"'’>1+1.3  <♦1+1.3  ■ ♦l.l’  - <“*>1+1.3  *■  <“>1+1.3  ■^<"‘>1.3]  ’ <'>•“> 
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■i+l.j 


-.j+lj  ’ 


(4)  f“ 


(C.20) 


^i+l.j  “ + -3"  (Ai^)  [-<““>1+1, j+i  [a^]  + 

■ ■'l.j+1  ’^i+lj+lj  , (C.21) 


(C.22) 


«i«,j  - ^[-  + (-'>-)i«.j]  * "j  {-  (j^) 

■>i.j+l^  " "±.j+l  '^i+lj+l  + ®l+l.j+l  - (a^] 

["i+i,j  ^n+i.j  - + ei+i,j  J + 

■ + (^)  { '"■’>l+l.J«  * - ^[-'i+i.j+1 

" 2 [-^+1,J  (ru),^j  + T,  J (ru)^^,  J j J (C.23) 


l+l,j  - '^o.^.  i+l,j  + 


(C.23) 


J 
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^1+1, j " 6 ^”'’^1+1, j ^i,j^  • 


®^l+l.j  “ -2  - "i"  (axJ  ^”“^1+1, j 


A 

‘'^1+1,  j = 2 


(C.25) 


(C.26) 


" -f ^ ["i+i,j  <n+i,j  - ♦!.])  + »i+i.jj- 

A •>/  (C.27) 

“ ■ 6 ^”'^^1+1, j+1  * (C.28) 


^1+l.j  “ 2 + 4"  [a^]  ^«“^i+l.j+l  • 

“6“  °^i+l,j+l  ■'■  ”i+l,j+l 

.»)}  ■ 


(C.29) 


A "1/2 

" 6 V®”  °^i+l,j+l  ■'■ 


- + ^i+1 


j+1  ^'*'i+l,j+l 


(C.30) 


A h * 


r ■ _ 


(C.31) 


(C.32) 
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where  the  various  parameters  in  the  above  equations  are  defined  by 


Wi+i,j  ^^®'‘^^i+i,j 

9 

(C.33) 

1 

fll  , fl]  1 

^i+l,j  “ 2Ax^ 

.Ui+i,j  Wi,jj  • 

(C.34) 

2 

2 2 

“ "e  i • 

■ “e  i+1  ■ “ l,j 

(C.35) 

1 

a - ^ 

*^1+1,  j Ax. 

X ^ 

W Ml)  1 ’ 

>■  -'i+l,j  ^^•'i,j-' 

(C.36) 

^o  i+1  *’  2 ^“i+1  “i^ 


(C.37) 


(C.38) 


'^i+l,j  " Ax^  i+l,j  ^i+l,j^  ■ Vl,j  ^"^i+l.j  ^i.j^  " '*'i,j^ 


- T 


i+l,j 


■ 

• ’ 

■ ''i,j 

% 

+ (s't) 


i.j 


i.jJ 


. (C.39) 


“i+l.j  Vl.j  ^''’ii-l.j  ■ "^ij^ 


(C.40) 


^i+l,j  “ ^^i+l,j  ^’^i.j  ■ ’*'i+l,j^  • (C.41) 

The  outer  region  must  be  considered  next.  Substitution  of  Equations 


(26),  (27),  and  (31)  into  Equation  (24)  yields  the  following  set  of  first- 
order,  nonlinear,  flnlte^lfference  expressions  for  the  outer  layer: 
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(rbT)^^l  - (rbT)j  - ^'^j+1^  ■•■  ~2  3x  j+1  “ % ^ 

-■^j  k ■*■  2^  12“  |-"j+l  k •*• 

3 ^ ^ j.  « 3 <■  2 2.  1 3 ,,  X fll+l  3 

3x  ^■'j+1^  2 3x  ^'*j+l  ■ '^e  ^ ' br^_^j^  3x  ^’‘’j+1^  [2  3x 


°'”j+i  - 'i+i  I?  ‘"j+i*  * 'j  I:  <'“>j  - <™)j  It  <’j>  - 1 h <“j^  - \'> 


1 3 X,  V fli  3 , 2 

bTJ^  [2  37 


/)  - abTj  - Tj  1^  (*j)J  } ’0  , (C.42) 


(Uj+i^  - - abT^^^ 


“j+l  - "l  - r*-  ''j+l  + ^j>  ■"  12-  'bJ^J-2^1?  <“l+l^  - “e^>  - ' 
-•J+1  ‘♦j+l>]  - [^1;  <“j^  - “e^>  -“'•^J  -^J  It  '♦j'  |-  0 


(C.43) 


Note  that  the  middle  equation  in  this  set  was  left  out  as  it  is  identical  to 
Equation  (C.2)  for  the  inner  layer.  The  same  procedure  is  applied  now  to  the 
above  equation  set  for  the  outer  region  as  was  applied  to  the  inner  set  of 
equations.  That  is,  they  are  approximated  in  the  center  of  the  box  using 
central  differences  and  then  are  linearized  using  Equations  (C.6)-(C.9).  The 
final  result  is  the  coupled  linear  system  for  the  corrections 
6u^^j^  j ’ ^^i+l  represented  by  Equations  {32)-(34).  The  form  is 

identical  to  that  of  the  inner  layer,  but  with  some  modifications  to  the 


coefficients  as  shown  below: 


A N ^1  r 1 

^^l+l,j  “ Ax^  ^"^l-H.j  6 [ax 


xj  °o  1+1, j |^“^'*'i+l,j  " '*'i,j^  AxJ 


^"l+l,j  "i.j^  ^o  i+1 


..1  ■ 


(C.44) 
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h ^ r 

..j  “ ■""f  (i;)  i"i.j  "i+i.j ' ^“Vi+i.j 


"■  ®o  i+i.j  ^’"i+i.j  - 


^i+i.j  “ -2  ^ ('^i+i,j  - + -i 


(C.45) 


(^)  |-2(ru),_ 


^0  i+i,J  ^’*'1+1, j " ’*'i,j^  ["^^Vi+l.j  ■ &x^ 


('l'i+1  j > » (C.46) 


■i+l,j  Ax, 


AxJ  ^^i+l.j+1  ^ 


^■^i+lj+l  ■*■  ■'i.J+l^  ~b~  [axJ  ®0  i+l,j+l  [^’*'1+1, j+1  " 

•1  ■*■  " ^0  i+l,j+l]  ’ (C.47) 


i+l,j  ” 3 


(ax^,  . 


j+1  ’^i+l.j+l  ■*■  ^“Vi+l,j+l 


^ = 2(rb),  ^ +.~^  ■*■ 

i+l,j  i+l,j+l  Ax^  ’^i+l,jri  i,j+l  o 


(C.48) 


“®0  i+l,j+l 


.j+1  “ ’’'i.j+l^  "^^Vi+l.j+l  ■ Ax^ 


^Vi.j+i  ■ '''i.j+iMr 


(C.49) 


February  2,  1979 
JMC : cac 


2 r 1 ' 

0 i+l,j+l  ■*■  ^i+l,j+l  i+l,j+l  “o  i+l,j  ■*■  ®i+i,j 


■^“o  ''0*1+1, j+1  ^*1+1, j+1  ■ '^l.j+l*  ■*■ 


*0  1+1  <"  1+1, 1+1  * ’'l+l,j+l'  [ ■'l+l,l+l  '"**1,1+1  * ■'1,1+1  <"**  1+1, l+l] 

" <<*0  'O*  1+1,1  ■ "l,!'  ■ "o  1+1  '“'l+l,!  ♦ '1+1,1*  * H"' 1+1,1 


'"*>1,1  'l,l  <"**1+1,1 


(C.50) 


6 [ax^J  ‘'o  i+l,j 


(C.51) 


(C.52) 


i+l,j 


--re 


* 

1+1, [a^]  ^’*'1+1, j ■ • 


(C.53) 


1+1, j 6 


(ax^)  1+1, 


j+1  ^^i+l,j+l  ■'■  "^l.j+l^  ’ (C.54) 


4 = 2 + -4“  (a„Bu) 


1+1, j 


3 1+1, j+1  * 


(C.55) 


•l+l,j 


1+1, j+1  {^*’^0^1+1, j+1  ■*■  (ax^  ^’''i+l.j+l  * ’‘'i,j+l^]  ’ 


(C.56) 


-42- 


February  2,  1979 
JMC : cac 


+ ^ 


i.j 


(C.57) 


where  6 


i+l,j’  ^0  i+1’  '^i,j  defined  as  before  and 


“o  I+l.J  ■ + ’i,j)  (l-irt.j  - . 


(C.58) 


0 i+l,j  “ 2 


Wi+i.j  * Wi.J  • 


"'O  m,j  - Si+I.j  (“i+i,/  + - A„  ^ + (b,)^_ J 

■ (s^)  “o  1+1.  j 


(C.59) 


The  remaining  coefficients  are  the  same  as  those  of  the  inner  layer. 


I 
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Appendix  D 

Block  Trldlagonal  Matrix  Solution 

This  appendix  describes  in  detail  the  solution  of  the  block  tridiagonal 
matrix  equation  given  by  Equation  (50) . 

Since  this  equation  is  block  tridiagonal,  L-U  factorization  can  be  used 
to  solve  it.  Following  Hoffman  [6], 

«=  LU  ; L = [Bj,  I,  0]  , U = [o,  a^,  (D.l) 

where  I is  the  unity  3x3  matrix  and  B,  a,  and  y are  3x3  matrices.  Equation 
(50)  then  can  be  expressed  as  the  two  systems 


LY  = R 


'I  ^ 

^ - 


UZ  = Y 


The  L-U  factorization  requires  that 


“l  = ^1 


and 


“ ®j  ♦ 2 < j < N+1 


Yj  = Cj  , 1 1 j < N , 


Yl  = . 


^j-1  ’ 2 < j < N+1 


“n+1  ^N+1  “ ^N+l  ’ 


■ ’'i  ■ ’'j  i ^ 1 


(D.2) 

(D.3) 

(D.4) 

(D.5) 

(D.6) 

(D.7) 

(D.8) 

(D.9) 

(D.IO) 

(D.ll) 
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Because  the  matrices  of  Equations  (42)-(46)  are  only  3x3,  the  linear 

systems  (D.4)-(D.ll)  can  be  Inverted  beforehand.  Thus  It  Is  not  necessary  to  use 

a generalized  matrix  Inversion  routine.  Equation  (D.4)  becomes 


10  0 
0 10 
A3,  ^ 


(D.12) 


L^l  “^1  -IJ 

For  j > 2 the  bottom  row  of  B is  zero;  thus,  from  Equation  (D.6)  the  bottom 
row  of  Bj  will  also  be  zero.  Then,  from  Equation  (D.5) 

“ll  “l2  “l3  ^j-1  $j-l  ^j-1 

“21  “22  “23  - %-!  ^j-1  ^j-1 

“31  “32  “33  . 


"11  "12 

®21  ^22  ®23 

0 0 

which  yields  the  relations 


®12  ®13 


0 0 
^1-1 


^“ll^j 

(“l2)j 

^“l3^j 

^“21^j 

^“22^j 

^“23^j 

(“3l)j 

(“32)j 

(“33)^ 


A , N 

°^j-l  ^ ^^13^j  °^j-l 

^^5-1  ' "Vl 

A , / A 

^j-1  ' ^®13^j  ^j-1 

^j-1  " ^®23^j  ^j-1 

^^j-1  " ^®23^j  ®^j-l 

^2j_i  - <^23^3  ^j-1 
A3j  . 

A 

■ 

“j  • 


(D.13) 


^ ^ * 
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I 


For  j=N+l,  Equation  (D.5)  becomes 


“11 

“12 

“13 

A 

A 

aN 

Aj, 

AN 

“21 

“22 

“23 

D2 

N 

F2 

N 

“31 

“32 

“33 

N+1 

0 

1 

0 

which  yields 


where  the  first  two  rows  are  the  same  as  Equations  (D.13)  with  j “ N+1  and 
j-l=N.  Only  the  last  row  changes  and  is  given  by  Eqxiation  (D.15).  ^ 

Two  cases  must  be  considered  for  Equation  (D.6),  namely,  j»2  and 
2 < j < N.  For  j=2. 


! ®11  ®12  ®13 

10  0 

^ A A 

Al-  Bl,  Cl, 

A \ A ^ 

j ®21  ®22  ^23 

jo  0 0 

2 

0 10 

II 

A2i  C2^ 

0 0 0 

(D.16) 


which  yields  the  solutions 

^2 


\ 


(612)2  " 


(6^3)2  - Ci,/A 


l'“2 


(623)2  " ^^1  ^1  ■ “1  ^1^/^ 
(622)2  ” ^3)/^ 


(D.17) 


(823)2  - C2^/A2 


f: 


- 

— 

p 

- 

1 

^11 

^12 

®13 

0 

0 

0 

^21 

CM 

CM 

OQ. 

623 

0 

0 

0 

(D.14) 

0 

0 

0 

N+1 

^N 

^N 

A 

■^N 

9 1 

i 

j 

1 

( 

“11 

“12 

“13 

N+1  ” 

®21 

®22 

“23 

1 

0 

1 

0 

N+1 

9 

(D.lb)  1 

:*• 


it 


liS' 
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i 

I 

1 

1 

3 


where 


A 

42  = C3^ 


(D.18) 


For  j > 2,  Equation  (D.6)  becomes 


®11  ^12  ^13 

^21  ®22  ®23 

“11  “12  “l3 
“21  “22  “23 

8B 

(i.  . M A 

A J-1  A 3-1  A j-1 

^^j-1  ®^j-l  ^^j-1 

, (0.19) 

0 

0 

0 J 

j 

O31  O32  033 

j-1 

0 0 0 

which  gives 
(6 


*23  “31^j-l 


022  “33^j-l  “21  “32^j-l  “2 

- (Cl  022  °31^j-l  ■ “23  “32^j-l  “ ®21  “33^j-lJ 

- r 

- (B1  Oj^  “i3^j-l  " (Ci  0^2  “ll^j-i  " (Al  0^3  “i2^j-l 

^®13^j  “ “11  “22^j-l  “21  “l3^j-l  ®12  “23^j-l 

^J-1  ■ “12  “21^j-l]  ‘ 

2i)j  - 022  “33^j-l  ■*■  “21  “32^j-l  ■*■  “23  “31^j-l 


KBl  <*11  a33)j_^  + (Cl  032  o^3)j_3  + (^1  03^  o^2>j-l 


/A. 


(A  O22  ®i3)j_i  - (B1 


>(0.20) 


(C2  O22  “3i)j_i  ■ (A2  O23  “32)j_i  “ (®2  02^^  “33^j_l 


/A 


j 


1(B2  O33  “33)j_i  + (A2  O32  033) j_2  + (C2  O33  032^3-! 


^®22^j  ” 

(B2  O33  033)j_3  - (C2  O32  “33  “i2^j-i 

(B23^j  ’ «11  “22^j-l  ■'■  “21  “l3^j-l  “12  “23^j-l 

* J 


(C  O22  »i3)j_i  - 0*33  “23^J-1  ■ “12  “21^J-1 
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where 


Aj  = 022  “33^j-l  ^“21  “32  “l3^j-l  ^“31  “l2  “23^j-l 


(031  O22  “ (“ii  “32  “23^j-l  ■ ^“21  “12  “3J^j-l  * 

Next,  solve  Equations  (D.8)  and  (D.9)  for  . Now 


(D.22) 


Hence,  Equation  (D.8)  yields  the  three  relations 


(Yl)^  = 0 , 


(Y2)  - 0^  , 

(Y3)^  = S3^  . 


Equation  (D.9)  reads 


^j-i  ^11  ^12  ^12 

5^j-l  ■ ^21  ^22  ®23  i 

S3^  0 0 0 j Y3j^_^ 

_J  L_  I—  — ' 


Upon  expansion,  for  2 < j < N+1,  the  solutions  are 


(D.23) 


(D.24) 


(Yl)j  = Slj_^  - (8^^)j  (Yl)^_^  - (Bj2^j  ■ ^®13>j  ’ 

A 

(Y2)^  = S2j_^  - (62^^  - (623)^  (Y3)j_^  , \(D.25) 

(V3)j  - , 


except  that  at  j“N+l, 


(V3)^j  • 0 . 


(D.26) 


\ 

* 


) 
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which  follows  from  Equation  (49).  Now  solve  Equation  (D.IO)  for  Z, 
using  Equation  (D.15): 


N+1 


“11 

®21 

“12  “l3 

®22  “23 

Sf 

Su 

Y1 

Y2 

1 

. ' ! 

(D.27)  i 

0 

M 

0 

1 

N+1 

St 

N+1 

0 

i 

N+1  ♦ i 

which  upon  expansion  and  solving  for  the  corrections  yields 


N+1 


P“l3Wl  ^^Wl  " ^“23^N+1  ^^^WiK®N+1 


= 0 


(D.28) 


(fit) 


W-1 


= ( 


^“2lWl  ^^^^N+1  ■ ^“ll^N+1  ^^^Wl  K ®N+1 


where 


■^N+l 


^“l3  “21W1  “ ^“11  “23^+1 


(D.29) 


Next,  solve  for  from  Equation  (D.ll)  for  2 < j < N. 

0 


“11 

“12 

“l3 

6i|i 

Y1 

“21 

“22 

“23 

6u 

» • 

Y2 

- 

“31 

O32 

“33 

j 

6t 

j 

Y3 
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which  yields  the  three  equations  for  the  correction  vectors 
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where 


r3j  = (Y3)j  - Aj  » 


(D.32) 


and  is  Identical  to  defined  by  Equation  (D.21),  or. 
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Note;  Equations  (D.31)  hold  for  2 < j < N.  For  j=l,  use  Equation  (D.12)  to  get 


which  yields  the  solutions 
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where  r3  is  given  by  Equation  (D.32)  and  ^2  Equation  (D.18). 
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